home *** CD-ROM | disk | FTP | other *** search
/ The Fatted Calf / The Fatted Calf.iso / Applications / Astronomy / Moon / Source / astro.c < prev    next >
Text File  |  1993-01-19  |  10KB  |  344 lines
  1. /* astro.c
  2.  * Part of the Moon application for the NeXT computer.
  3.  * Authors:  John Walker of Autodesk
  4.  *           Geoffrey S. Knauth (NeXT port)
  5.  * Date:     January 4, 1992
  6.  *
  7.  * This code was placed in the public domain by John Walker.
  8.  * phasehunt has been rewritten.  Define USE_GSK_PHASEHUNT to use it.
  9.  */
  10.  
  11. #import "all.h"
  12.  
  13. #define TRUE    1
  14. #define FALSE    0
  15.  
  16. /*  MEANPHASE  --  Calculates mean phase of the Moon for a given
  17.            base date and desired phase:
  18.                0.0   New Moon
  19.                0.25  First quarter
  20.                0.5   Full moon
  21.                0.75  Last quarter
  22.             Beware!!!  This routine returns meaningless
  23.                     results for any other phase arguments.  Don't
  24.             attempt to generalise it without understanding
  25.             that the motion of the moon is far more complicated
  26.             that this calculation reveals. */
  27.  
  28. static double meanphase(sdate, phase, usek)
  29. double sdate, phase;
  30. double *usek;
  31. {
  32.     int yy, mm, dd;
  33.     double k, t, t2, t3, nt1;
  34.  
  35.     jyear(sdate, &yy, &mm, &dd);
  36.  
  37.     k = (yy + ((mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685;
  38.  
  39.     /* Time in Julian centuries from 1900 January 0.5 */
  40.     t = (sdate - 2415020.0) / 36525;
  41.     t2 = t * t;           /* Square for frequent use */
  42.     t3 = t2 * t;           /* Cube for frequent use */
  43.  
  44.     *usek = k = floor(k) + phase;
  45.     nt1 = 2415020.75933 + synmonth * k
  46.           + 0.0001178 * t2
  47.           - 0.000000155 * t3
  48.           + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
  49.  
  50.     return nt1;
  51. }
  52.  
  53. /*  TRUEPHASE  --  Given a K value used to determine the
  54.            mean phase of the new moon, and a phase
  55.            selector (0.0, 0.25, 0.5, 0.75), obtain
  56.            the true, corrected phase time.  */
  57.  
  58. static double truephase(k, phase)
  59. double k, phase;
  60. {
  61.     double t, t2, t3, pt, m, mprime, f;
  62.     int apcor = FALSE;
  63.  
  64.     k += phase;           /* Add phase to new moon time */
  65.     t = k / 1236.85;       /* Time in Julian centuries from
  66.                       1900 January 0.5 */
  67.     t2 = t * t;           /* Square for frequent use */
  68.     t3 = t2 * t;           /* Cube for frequent use */
  69.     pt = 2415020.75933       /* Mean time of phase */
  70.          + synmonth * k
  71.          + 0.0001178 * t2
  72.          - 0.000000155 * t3
  73.          + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
  74.  
  75.         m = 359.2242               /* Sun's mean anomaly */
  76.         + 29.10535608 * k
  77.         - 0.0000333 * t2
  78.         - 0.00000347 * t3;
  79.         mprime = 306.0253          /* Moon's mean anomaly */
  80.         + 385.81691806 * k
  81.         + 0.0107306 * t2
  82.         + 0.00001236 * t3;
  83.         f = 21.2964                /* Moon's argument of latitude */
  84.         + 390.67050646 * k
  85.         - 0.0016528 * t2
  86.         - 0.00000239 * t3;
  87.     if ((phase < 0.01) || (abs(phase - 0.5) < 0.01)) {
  88.  
  89.        /* Corrections for New and Full Moon */
  90.  
  91.        pt +=     (0.1734 - 0.000393 * t) * dsin(m)
  92.             + 0.0021 * dsin(2 * m)
  93.             - 0.4068 * dsin(mprime)
  94.             + 0.0161 * dsin(2 * mprime)
  95.             - 0.0004 * dsin(3 * mprime)
  96.             + 0.0104 * dsin(2 * f)
  97.             - 0.0051 * dsin(m + mprime)
  98.             - 0.0074 * dsin(m - mprime)
  99.             + 0.0004 * dsin(2 * f + m)
  100.             - 0.0004 * dsin(2 * f - m)
  101.             - 0.0006 * dsin(2 * f + mprime)
  102.             + 0.0010 * dsin(2 * f - mprime)
  103.             + 0.0005 * dsin(m + 2 * mprime);
  104.        apcor = TRUE;
  105.     } else if ((abs(phase - 0.25) < 0.01 || (abs(phase - 0.75) < 0.01))) {
  106.        pt +=     (0.1721 - 0.0004 * t) * dsin(m)
  107.             + 0.0021 * dsin(2 * m)
  108.             - 0.6280 * dsin(mprime)
  109.             + 0.0089 * dsin(2 * mprime)
  110.             - 0.0004 * dsin(3 * mprime)
  111.             + 0.0079 * dsin(2 * f)
  112.             - 0.0119 * dsin(m + mprime)
  113.             - 0.0047 * dsin(m - mprime)
  114.             + 0.0003 * dsin(2 * f + m)
  115.             - 0.0004 * dsin(2 * f - m)
  116.             - 0.0006 * dsin(2 * f + mprime)
  117.             + 0.0021 * dsin(2 * f - mprime)
  118.             + 0.0003 * dsin(m + 2 * mprime)
  119.             + 0.0004 * dsin(m - 2 * mprime)
  120.             - 0.0003 * dsin(2 * m + mprime);
  121.        if (phase < 0.5)
  122.           /* First quarter correction */
  123.           pt += 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime);
  124.        else
  125.           /* Last quarter correction */
  126.           pt += -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime);
  127.        apcor = TRUE;
  128.     }
  129.     if (!apcor) {
  130.            fprintf(stderr, "TRUEPHASE called with invalid phase selector.\n");
  131.        abort();
  132.     }
  133.     return pt;
  134. }
  135.  
  136. /*  PHASEHUNT  --  Find time of phases of the moon which surround
  137.            the current date.  Five phases are found, starting
  138.            and ending with the new moons which bound the
  139.            current lunation.  */
  140.  
  141. #define USE_GSK_PHASEHUNT
  142. #ifdef    USE_GSK_PHASEHUNT
  143.  
  144. void phasehunt(sdate, phases)
  145. double sdate;
  146. double phases[5];
  147. {
  148.     double adate, nt, k, k1, k2;
  149.     double nextNewMoon = 0.0;
  150.  
  151.   /* Go back two months and about a second. */
  152.     adate = sdate - (synmonth * 2.00001);
  153.  
  154.   /* Awful meanphase has caused me major headaches.
  155.    * It skips lunations, and provides only the roughest estimate.
  156.    * I will only call it to get the initial seed value for k.
  157.    */
  158.     nt = meanphase(adate, 0.0, &k);
  159.  
  160.     for (k2 = k; sdate > nextNewMoon; k += 1.0)
  161.     nextNewMoon = truephase(k2 = k, 0.0);
  162.     k1 = k2 - 1.0;
  163.  
  164.     phases[0] = truephase(k1, 0.0);    /* new moon */
  165.     phases[1] = truephase(k1, 0.25);    /* first quarter */
  166.     phases[2] = truephase(k1, 0.5);    /* full moon */
  167.     phases[3] = truephase(k1, 0.75);    /* last quarter */
  168.     phases[4] = nextNewMoon;        /* next new moon */
  169. }
  170.  
  171. #else    !USE_GSK_PHASEHUNT
  172.  
  173. static void phasehunt(sdate, phases)
  174. double sdate;
  175. double phases[5];
  176. {
  177.     double adate, k1, k2, nt1, nt2;
  178.  
  179.     adate = sdate - 45;
  180.     nt1 = meanphase(adate, 0.0, &k1);
  181.     while (TRUE) {
  182.        adate += synmonth;
  183.        nt2 = meanphase(adate, 0.0, &k2);
  184.        if (nt1 <= sdate && nt2 > sdate)
  185.           break;
  186.        nt1 = nt2;
  187.        k1 = k2;
  188.     }
  189.     phases[0] = truephase(k1, 0.0);
  190.     phases[1] = truephase(k1, 0.25);
  191.     phases[2] = truephase(k1, 0.5);
  192.     phases[3] = truephase(k1, 0.75);
  193.     phases[4] = truephase(k2, 0.0);
  194. }
  195.  
  196. #endif    !USE_GSK_PHASEHUNT
  197.  
  198. /*  KEPLER  --    Solve the equation of Kepler.  */
  199.  
  200. static double kepler(m, ecc)
  201. double m, ecc;
  202. {
  203.     double e, delta;
  204. #define EPSILON 1E-6
  205.  
  206.     e = m = torad(m);
  207.     do {
  208.        delta = e - ecc * sin(e) - m;
  209.        e -= delta / (1 - ecc * cos(e));
  210.     } while (abs(delta) > EPSILON);
  211.     return e;
  212. }
  213.  
  214. /*  PHASE  --  Calculate phase of moon as a fraction:
  215.  
  216.     The argument is the time for which the phase is requested,
  217.     expressed as a Julian date and fraction.  Returns the terminator
  218.     phase angle as a percentage of a full circle (i.e., 0 to 1),
  219.     and stores into pointer arguments the illuminated fraction of
  220.         the Moon's disc, the Moon's age in days and fraction, the
  221.     distance of the Moon from the centre of the Earth, and the
  222.     angular diameter subtended by the Moon as seen by an observer
  223.     at the centre of the Earth.
  224.  
  225. */
  226.  
  227. double jdtophase(pdate, pphase, mage, dist, angdia, sudist, suangdia)
  228. double pdate;
  229. double *pphase;            /* Illuminated fraction */
  230. double *mage;               /* Age of moon in days */
  231. double *dist;               /* Distance in kilometres */
  232. double *angdia;            /* Angular diameter in degrees */
  233. double *sudist;            /* Distance to Sun */
  234. double *suangdia;                  /* Sun's angular diameter */
  235. {
  236.  
  237.     double Day, N, M, Ec, Lambdasun, ml, MM, MN, Ev, Ae, A3, MmP,
  238.            mEc, A4, lP, V, lPP, NP, y, x, Lambdamoon, BetaM,
  239.            MoonAge, MoonPhase,
  240.            MoonDist, MoonDFrac, MoonAng, MoonPar,
  241.            F, SunDist, SunAng;
  242.  
  243.         /* Calculation of the Sun's position */
  244.  
  245.     Day = pdate - epoch;        /* Date within epoch */
  246.     N = fixangle((360 / 365.2422) * Day); /* Mean anomaly of the Sun */
  247.     M = fixangle(N + elonge - elongp);    /* Convert from perigee
  248.                        co-ordinates to epoch 1980.0 */
  249.     Ec = kepler(M, eccent);     /* Solve equation of Kepler */
  250.     Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
  251.     Ec = 2 * todeg(atan(Ec));   /* True anomaly */
  252.         Lambdasun = fixangle(Ec + elongp);  /* Sun's geocentric ecliptic
  253.                            longitude */
  254.     /* Orbital distance factor */
  255.     F = ((1 + eccent * cos(torad(Ec))) / (1 - eccent * eccent));
  256.     SunDist = sunsmax / F;        /* Distance to Sun in km */
  257.         SunAng = F * sunangsiz;     /* Sun's angular size in degrees */
  258.  
  259.  
  260.         /* Calculation of the Moon's position */
  261.  
  262.         /* Moon's mean longitude */
  263.     ml = fixangle(13.1763966 * Day + mmlong);
  264.  
  265.         /* Moon's mean anomaly */
  266.     MM = fixangle(ml - 0.1114041 * Day - mmlongp);
  267.  
  268.         /* Moon's ascending node mean longitude */
  269.     MN = fixangle(mlnode - 0.0529539 * Day);
  270.  
  271.     /* Evection */
  272.     Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
  273.  
  274.     /* Annual equation */
  275.     Ae = 0.1858 * sin(torad(M));
  276.  
  277.     /* Correction term */
  278.     A3 = 0.37 * sin(torad(M));
  279.  
  280.     /* Corrected anomaly */
  281.     MmP = MM + Ev - Ae - A3;
  282.  
  283.     /* Correction for the equation of the centre */
  284.     mEc = 6.2886 * sin(torad(MmP));
  285.  
  286.     /* Another correction term */
  287.     A4 = 0.214 * sin(torad(2 * MmP));
  288.  
  289.     /* Corrected longitude */
  290.     lP = ml + Ev + mEc - Ae + A4;
  291.  
  292.     /* Variation */
  293.     V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
  294.  
  295.     /* True longitude */
  296.     lPP = lP + V;
  297.  
  298.     /* Corrected longitude of the node */
  299.     NP = MN - 0.16 * sin(torad(M));
  300.  
  301.     /* Y inclination coordinate */
  302.     y = sin(torad(lPP - NP)) * cos(torad(minc));
  303.  
  304.     /* X inclination coordinate */
  305.     x = cos(torad(lPP - NP));
  306.  
  307.     /* Ecliptic longitude */
  308.     Lambdamoon = todeg(atan2(y, x));
  309.     Lambdamoon += NP;
  310.  
  311.     /* Ecliptic latitude */
  312.     BetaM = todeg(asin(sin(torad(lPP - NP)) * sin(torad(minc))));
  313.  
  314.     /* Calculation of the phase of the Moon */
  315.  
  316.     /* Age of the Moon in degrees */
  317.     MoonAge = lPP - Lambdasun;
  318.  
  319.     /* Phase of the Moon */
  320.     MoonPhase = (1 - cos(torad(MoonAge))) / 2;
  321.  
  322.     /* Calculate distance of moon from the centre of the Earth */
  323.  
  324.     MoonDist = (msmax * (1 - mecc * mecc)) /
  325.        (1 + mecc * cos(torad(MmP + mEc)));
  326.  
  327.         /* Calculate Moon's angular diameter */
  328.  
  329.     MoonDFrac = MoonDist / msmax;
  330.     MoonAng = mangsiz / MoonDFrac;
  331.  
  332.         /* Calculate Moon's parallax */
  333.  
  334.     MoonPar = mparallax / MoonDFrac;
  335.  
  336.     *pphase = MoonPhase;
  337.     *mage = synmonth * (fixangle(MoonAge) / 360.0);
  338.     *dist = MoonDist;
  339.     *angdia = MoonAng;
  340.     *sudist = SunDist;
  341.     *suangdia = SunAng;
  342.     return fixangle(MoonAge) / 360.0;
  343. }
  344.